(* Mathematica package *)
Options[GenerateInterfacialModel] = {"ReplaceRules" -> {}, "ConstraintRules" -> {}};

GenerateInterfacialModel[
  dataset_Integer, 
  InputTransferCArgs : {{{(*Ep*)_,(*Es*)_}, (*Omega*)_, (*Theta0*)_List, (*n*)_List, (*d*)_List, (*nbar*)_List}..},
  UnitChi_,
  OutputTransferCArgs:{(*OmegaWM*)_, (*ThetaPM*)_List, (*nlist*)_List, (*d*)_List, (*nbarINT*)_List},
  {numD_, numV_, numR_},  
  Parameters : {M_, {omega0_, A_, gamma_}, ANR_}, OptionsPattern[]] :=
 
 Module[{positions = Position[UnitChi, Except[0, _Integer]], 
   elements, ChiFunc, Chi}, 
  elements = 
   StringJoin /@ (positions /. {1 -> "x", 2 -> "y", 3 -> "z"});

  
  
  ChiFunc = (Sqrt[Abs[M[dataset]]] * Apply[Plus, 
  	LorentzianOSC[OutputTransferCArgs[[1]], Table[head[ele, dataset,v,r],{head,{omega0,A,gamma}},{ele, elements},{v,numV},{r,numR}]],{2}] +
  	Table[ANR[ele,dataset,v],{ele,elements},{v,numV}]) /. OptionValue["ReplaceRules"] /. OptionValue["ConstraintRules"];
  
  Table[Chi[v] = ReplacePart[UnitChi, Table[
  positions[[i]] -> ChiFunc[[i, v]], {i, Length @ positions}]], {v, numV}];
  
  Sum[KroneckerDelta[d-dataset] * Table[OutputTransferC[#, Sequence @@ (MapAt[Sequence@@{Part[#, v],v} &, #, -1] & @ OutputTransferCArgs)]& @
     					(NonlinearPolarization[Chi[v], 
     										   InputTransferC @@@ (MapAt[Sequence@@{Part[#, v],v} &, #, -1] & /@ InputTransferCArgs)]), {v, numV}],{d,numD}]
  ]

GenerateInterfacialModel[
  InputTransferCArgs : {{{(*Ep*)_,(*Es*)_}, (*Omega*)_, (*Theta0*)_List, (*n*)_List, (*d*)_List, (*nbar*)_List}..},
  ChiV_,
  OutputTransferCArgs:{(*OmegaWM*)_, (*ThetaPM*)_List, (*nlist*)_List, (*d*)_List, (*nbarINT*)_List}] :=
	Table[
		With[{pNL = NonlinearPolarization[ChiV[[v]], InputTransferC @@@ (MapAt[Sequence@@{Part[#, v],v} &, #, -1] & /@ InputTransferCArgs)]},

		OutputTransferC[
  						pNL,
						Sequence @@ (MapAt[Sequence@@{Part[#, v],v++} &, OutputTransferCArgs, -1])]], {v, Length@ChiV}]


(*ChiFunc = (Sqrt[Abs[M[dataset]]] * Apply[Plus, 
  	LorentzianOSC[OutputTransferCArgs[[1]], Table[head[ele, dataset,v,r],{head,{omega0,A,gamma}},{ele, elements},{v,numV},{r,numR}]],{2}] +
  	Table[ANR[ele,dataset,v],{ele,elements},{v,numV}]) /. OptionValue["ReplaceRules"] /. OptionValue["ConstraintRules"]*)